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Our interest in averaging techniques steins from a requirement to conduct large-scale mission 
planning exercises on a system of computers that gets overloaded whenever si*ch a study 
arises. These studies have been conducted at Aerospace using a Cowell propagator simply 
because of the necessity for accuracy in a high drag environment. However, due to limited 
computer resources^ a project has been initiated to develop a better method, which will 
involve using an existing algorithm or developing a new one. 

This paper presents two topics, tne first of which concerns a feasibility study that was con- 
ducted using two programs, MAESTRO and PECOS, to, see if averaging could be accurate. 
The second topic concerns some continuing investigations and some new ideas, generally 
untested. It is assumed that the idea of averaging with respect to trajectory propagation is 
familiar to everyone so it wiL not be reviewed here. 

The Maestro program, developed by Chauncey Uphoff and Dave Lutzky of Analytical 
Mccimnics Associates, was obtained by Aerospace for a feasibility test. It was initially 
designed as an interplanetary mission analysis tool for the Radio Astronomy Explorer-2 
(RAE-2) lunar orbiter. Several choices of variables and propagating techniques are available 
in M\ESTRO. For instance, a Cowell integrator, an Encke integrator, and several variation- 
of-parametor techniques, deluding the numerical averaging of the Gaussian variation-of- 
panameter ( VOP) equations, are available. 

The forces in the program are given in the radial, circumferential, and polar directions. 
MAESTRO has a fairly simple exponential atmosphere and has a precision averaging startup, 
wiiich means that osculating initial conditions in several choices of reference frames can be 
input and the program will automatically provide the initial mean state. 

The Aerospace version of MAESTRO has been merged with a program called PECOS. 

PECOS, itself, is a test bed for orbit planning simulations. Tlic interplanetary calculations 
were stripped out of MAESTRO and the force evaluations from PECOS are now used by 
both MAESTRO and PECOS. 

The output can be plotted by an auxiliary pr a i available at Aerospace. In particular, 
overlays of PECOS p.^-ecision and MAESTRO averaging output can easily be accomplished for 
comparison purposes. Two versions of mean-to-osculatinir * ar'* u ed. 
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On output, the Kozai form is used, and within the averaging quar'rature, the Iszak form is 
used. There is no particular reason for this order except for the program’s historical develop- 
ment. Fortunately, Iszak is very good for use under the quadrature because it is nonsinguiar 
at low eccentricity. 

With present program limitations, only certain parameter sets and integration types from 
MAESTRO are available. For instance, the Cowell integration is not available. Only two out 
of the original eight techniques are being used. They include the integration of a parameter 
set which is nonsingular at low eccentricity. This set is similar to. but not exactly the same 
as, the equinoctial elements. The averaging technique uses the same olcaients. 

The initial conditions for the test case that were used for the feasibility study were chosen 
from ty^acal orbits of interest to Aerospace. They were initially compiled for Terry Harter 
to run at GSFC on a test using the Goddard Trajectory Determination System (GTDS). At 
Aerospace they were run by Stan Navickas in the Mission Analysis Department with the 
MAESTRO/PECOS software. 

Table 1 shows that the four cases here have perigees ranging from 67 to 92 n.m. and apogees 
ranging from 1 56 to 223 n.m. The inclinations are all between 90^ and 1 1 0.5°. The life- 
times range from about 4.S days to 2 months; case A-1 has a 2nnonth lifetime and case 
B-l has a 4.S-day lifetime, and those are the two ca.ses that will be discussed here. 

Figures 1 . 2, and 3 are time histories of elements for case A-l , which had approximately a 
2-month lifetime. The computer runs were limited to 7 days because it is so expensive to 
run the Cowell integration in PECOS. This is the orbit which does the least of the four 
test cases. 


Table 1 

Initial Conditions for Four Test Cases 


Case 

Perigee 

Hp(n.m.) 

Apogee 

inciinauon (deg) 

C„ A/W(fr /lb) 

A1 

92 

156 

96.43 

0.0007295 

(0.002099) 

\2 

85 

I6I 

94.55 

0.0007295 

(0.002099) 

B1 

61 

207 

96.57 

0.0062981 

(0.001286) 

B2 

69 

223 

110 5 

0.0062981 

(0.001286) 
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The semimajor axis is plotted in figure I . The oscillating line is the output from the Cowell 
technique in PECOS, and the line down the center is the mean element output from 
MAESTRO. There has been no correction to put the short-period variations back on. 

Figure 2 is the eccentricity from the same trajectory. Again the oscillating curve is the out- 
put from the Cowell integration translated into Kepler elements. Figure 3 is the ai^ument 
of perigee. Figures 4, 5, and 6 are for case B-1 , which dies ir* 4.5 days. The dark l*ne down 
the center is the mean output from MAESTRO. It can be seen thai the program tracks the 
elements all the way to the end of the lifetime. 

Figure 6 is the argument of perigee, and, as expected, the envelope of the oscillations grows 
rapidly as it approaches the end of the lifetime and the eccentricity drops. The orbit be- 
comes circular, so the argument of perigee begins to osdliate more and more wildly. But 
the mean argument of perigee is still being tracked and follows the center of the oscillation 
envelope. 

Figures 7, 8, and 9 give the comparison between the Cowell integration and the corrected 
average propagation at about 4 days into the case having a 4.5-day lifetime. The orbit is 
generally tracked. The phasing and the magnitude have tracked quite well at the end of 4 
days. I am very hesitant to interpret the small errors there, because we have a problem with 
the interpolator in PECOS/MAESTRO right now. The step sizes in the integration for the 
averaging process are 90 minutes, but the output is roughly every 2 minutes. 

if the output of mean elements were plotted on a finer scale, it would be seen tnat the result 
has a small wiggle. This wiggle is probably due to the interpolator. There is some suspicion 
that the errors in the osdliating elements are due to that wiggle, but that is still a qualitative 
judgment. The iinportant thing is that the orbit is being tracked. 

We int^erpreted these plots as proof that the averaging technique does work. At this point, 
there had been no attempt whatsoever to make this program efficient. MAESTRO was sim- 
ply merged into a larger system using some of the logic from that larger system, so timinp 
from these runs were totally meaningless. Using the results proving the feasibility of the 
technique, we decided to go ahead and attempt to modify MAESTRO to make it a stand- 
alone version that would be efficient, using essentially the same logic that it has always had, 
but using a more accurate force model than it had when it was developed for GSFC. 

The modifications and developments which will go into this further study will be discussed 
now. 

The purpose of the current investigation is to increase the accuracy a~-i, if possible, to de- 
crease the computation time; this will of course involve a compromise between the two. 
CuTently, we are developing MAESTRO as a stand-alone version with modification for 
accuracy in a high drag environment, to run on the CDC 7600. The next step will be to 
replace the independent variable with the true anomaly and to remove the fast variable. 
Another concurrent activity, being handled by Dave Lutzky of Vector Sciences, is treating 
the short period effects within the averaging quadrature and for output with a Fourier ex- 
pansion from which the coefficients of the harmonics are computed automatically and 
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numericaily. This is work which was developed under funds from GSFC. The results will be 
inserted into MAESTRO to see if they will hdp to make the program more efficient and 
accurate. 

Now I would like to discuss some of the work I am doing. It is all formal at this time, none 
of it has been implemented in the program yet. 

One of the disadvant^es of using time as the independent parameter is that the averaging 
period is always difficult to interpret. For instance, it is possible to use the mean motion, 
which is determined by the initial oscillating elements or by the initial mean elements, to 
determine the period. But, however it is done, it will probably be difficult to be consistent 
and to initially choose the right time. 

Another problem is that the equations for the averaged elements have a term in their exact 
form which depends upon the time rate of the period. This is generally referred to as the 
Leibnitz term, and it depends upon the time rate of the period, the values of the true 
elements at both ends of the averaging period, and on the mean element at the midpoint of 
the averaging interval. 

Still another problem arises if constant time steps are chosen; For instance, if the integration 
step is chosen as the initial period of the orbit, then, farther on down the line, the time step 
will no longer correspond to the period, and there will be drift with respect to perigee or any 
geometric event in the orbit. If the true anomaly is chosen as the independent parameter, 
none of the above disadvantages apply. Another effect is that we get automatic regulariza- 
tion of the orbit sampling. Points are automatically clustered near perigee where most of the 
action occurs for the high drag situation. 

It is also possible to make repeated use of terms in the mean-to-osculating transfcH'mations, 
particularly if step size is chosen properly within the averaging quadrature. Also, the mean- 
to-osculating formulas appear in closed form, at least to first order in . This is due simply 
to the fact that the disturbing potential can be expressed in a finite number of terms, de- 
pending upon the true anomaly. 

In changing the propagation algorithm, we wanted to choose a state vector which is nonsin- 
gular for low eccentricity and inclination. For convenience, we chose to use the same state 
that MAESTRO now uses, except for the fast variable. If the independent variable is now 
chosen as the true anomaly, then the angle-time relationship must be tracked through an 
equation other than the one we had before. If the new fast variable is itself chosen as the 
time, we then have a new differential equation; 

dt - 

— = h/r + perturbations, 
df 

Because this leaves us with large oscillations, it is a poor choice for averaging. Suppose 
that, instead, the parameter Q (Stem, 1960) is chosen; 
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Q = M + cj + S2-nt 
n = (-2c)^ III 


where e is the orbital energ>’. Then h determines the period of the true anomaly. The Q 
will always be a slow variable, and, in fact, for the Kepler problem, it is a constant. 

The equations for the new system are given by 


wherj L is any element 


IE dE dt 
df ~dt df 


(E,t,0 


A 

df 


E=F„(E.t,0 


The exact averaged equations are 


E = 


- r 

J 

f-W 


EfOdf* 


df 



dt 


Fg(E,r,t)df'~ 


To make averaging efficient, the force evaluation is accomplished by approximating the 
elements in the int^rand to the mean elements, plus a correction, as in the standard second- 
order averaging technique. Also, since the time no longer appears explicitly in the state, it 
must be the result of a similar correcting transformation. This means that the mean-to- 
osculating transformation is required at every step, as it is in all second-order averaging 
techniques. 

Formally, the mean-to-osculating transformations are simple to derive. The equations 
needed are indicated below in abbreviated form: 


E* = E + 6Ep 


If it is assumed that Af is 2rr/2K, then values of the tngonometric functions of the harmon- 
ics of the true anomaly can be tabulated and never recalculated, because the same values 
will be repeated in true anomaly from revolution to revolution. This should save a signifi- 
cant amount of calculation and make for an efficient mean-to-osculating transformation. 


In the equations, 


N 


dE 

— = CONST + 
df 



(A cos nf + sin nO 

' n n 
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and 


N 

6Ep =^^(-A^ sin nf + cos nO/n, 

I 

and are derived from the Lagrange planetary equations with the J ^ perturbation only. 
The ( J are slowly varying functions of the state vector and supposedly would not 
need to be recalculated more than once per revolution, and probably not that often. This 
depends upon the strength of the perturbations. 

The evaluation of dQ/df is fairly obvious, and it is used here as an example; The forcing 
function, , is (M + w + S2 - n) + (n - n) - nt. Time appears explicitly, which may create 
some difficulties, at least formally, although it is doubtful that numerically this will give 
any trouble. The factor dt/df, which appears in the integrand, is given by 

dt h 

— = — + perturbations, 
df 

However, since all the parameters in the state are slow (there is no fast variable), probably 
the two-body rate, or just h/r^ , is good enough for the function dt/df. That will probably 
require some experimentation, but, since there is no fast variable, there is no purpose in 
carrying higher order terms in this factor. The term h is related to the time derivative of 
the eneigy, and that relates (at least in its dominant term) to the time rate of the semimajm' 
axis due to drag for the cases where drag is the most impcvtant perturbation aside from . 
The quantity (n - n) is fH-oportional to the disturbing potential and (M cu -t- A - n) is 
given by the Gaussian VOP equations. The MAESTRO propagator will be used as the basic 
software tool in the implementation of these techniques. 

Still to be considered are techniques for calculating an accurate argument of perigee at low 
eccentricity. There are some problems because there is a mix of the true mean anomaly 
and the mean argument of perigee in one of the transformations; some experimentation is 
still required. 

The exponential instability in Q, which is due to the presence of ht, must also be investigated. 
This is a problem which will most likely be more apparent on computers with short word 
lengths. With the CDC computers and their essentially inflnite word lengths, we have tried 
experimenting to see when we can make the equations blow up in a simple form, and we 
have been unable to And the instability numerically. That may give us trouble, though, 
and it will require more experimentation to see if there really is a problem. 
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Figure 5. Comparison of PECOS and (MAESTRO eccentricity plots for case B-1 
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Figure 6. Comparison of PECOS and MAESTRO argument of perigee plots for case B*l. 
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Figure 8. Comparison of Cowell integration and short period recovery for eccentricity. 
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